function [ A, b, z0 ] = makeLCP( obj )
    M = obj.dynamics.M;
    Minv = obj.dynamics.Minv;
    Gn = obj.dynamics.Gn;
    Gf = obj.dynamics.Gf;
    E  = obj.dynamics.E;
    U = obj.dynamics.U;

    NU = obj.dynamics.Vel;
    nb = size(M, 1) / 6;
    nc = length(obj.dynamics.PSI);
    nd = obj.dynamics.nd;
    if isfield(obj, 'constraints')
        % bilateral constraints
        nd = length(obj.constraints.pairs);
    end
    h = obj.dynamics.h;
    FX = obj.dynamics.Forces;
    PSI = obj.contacts.PSI;   % the depth is actually the penetration here 
%      figure
%      spy(Minv);
%      title('Minv');
%      figure
%      spy(Gn);
%      title('Gn')
    MinvGn = Minv * Gn;
    MinvGf = Minv * Gf;
    
    A = [ Gn'*MinvGn   Gn'*MinvGf  zeros(nc);
          Gf'*MinvGn   Gf'*MinvGf  E;
          U            -E'         zeros(nc)];
    ASYD = [ Gn'*MinvGn   Gn'*MinvGf  ;
          Gf'*MinvGn   Gf'*MinvGf  ];
 
 
    b = [ Gn'*(NU + M\FX*h) + PSI/h;    % FX*h could be replaced if we stored Pext instead of Fext
              Gf'*(NU + M\FX*h);
                zeros(nc,1) ];
    
    z0 = zeros(length(b),1); 
end

